

clear

rng(328120)

rawdata = csvread('../../01_data/02_basedata/estimationdata/All_T_J/All_T_J_choices_buff_0_0_lf_principal.csv',1,0);



IJVAmadev       = rawdata(:,4);
IJVAmean        = rawdata(:,6);

VA = IJVAmadev+IJVAmean;

clearvars rawdata;

load '../../01_data/05_temp/mixed_logit_principal_noVA'    

IJX = IJX(VA~=0,:);
IJjtidx = IJjtidx(VA~=0,:);
Xmat = Xmat(VA~=0,:);


error_draws=evrnd(0,1,size(IJX,1),1); 

IJjpref=Xmat*thetahat(1:end-1) + error_draws;

% rank positions
pref_temp = IJjpref;

rank_vec = zeros(size(pref_temp));


INC = accumarray(IJjtidx,ones(size(pref_temp)),[],@sum);

maxC = max(INC);

for cc=1:maxC
    maxcurrent = accumarray(IJjtidx,pref_temp,[],@max);
    maxcurrent = max(maxcurrent,-1000000+1);
    rank_vec(pref_temp==maxcurrent(IJjtidx)) = cc;
    pref_temp(pref_temp==maxcurrent(IJjtidx)) = -1000000;
end

rank_vec = rank_vec ./ INC(IJjtidx);

rank_vec = 100*(1-rank_vec);

% output

outputvec = VA(VA~=0);
rankvecoutput = rank_vec;

output_q =  quantile(outputvec,20);

output_v = zeros(size(outputvec));
for vv=1:20 
   output_v(output_v==0 & outputvec<=output_q(vv)) = vv; 
end
output_v(output_v==0 & outputvec>output_q(20)) = vv+1;

mean_rank_output = accumarray(output_v,rankvecoutput,[],@mean);
num_rank_output = accumarray(output_v,ones(size(rankvecoutput)),[],@sum);
output_ceil_short = accumarray(output_v,outputvec,[],@mean);

output_ceil_short(num_rank_output==0) = [];
mean_rank_output(num_rank_output==0) = [];
p10_rank_output = zeros(21,1);
p90_rank_output = zeros(21,1);
for vv=1:21
    p10_rank_output(vv) = prctile(rank_vec(output_v==vv),10);
    p90_rank_output(vv) = prctile(rank_vec(output_v==vv),90);
end
p10_rank_output(num_rank_output==0) = [];
p90_rank_output(num_rank_output==0) = [];

clf

figure(1)
scatter(output_ceil_short(2:end),p10_rank_output(2:end),'x');
hold on
scatter(output_ceil_short(2:end),mean_rank_output(2:end),'o');
hold on
scatter(output_ceil_short(2:end),p90_rank_output(2:end),'+');
xlabel('total value added')
ylabel('preference percentile')
ylim([0 100])
saveas(gcf,'../../03_output/figures/principal_rank_output_noVA.png')


